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The problem of how to compute accurately and efficiently the sign of the overlap between two 
general HFB wave functions is addressed. The results obtained can easily be extrapolated to the 
evaluation of the sign of the trace of a density operator exponential of one body operators. 



I. INTRODUCTION 

Beyond mean field calculations are becoming very pop- 
ular [l|] as they allow a fairly good description of many 
nuclear state properties of both the ground state and sev- 
eral kinds of excited states all over the Nuclide chart. In 
these calculations, overlaps of Hartree- Fock- Bogoliubov 
(HFB) wave functions have to be computed. Standard 
formulas [2] involve the square root of a determinant leav- 
ing the sign of the overlap undefined. However, when the 
HFB states preserve some kind of discrete symmetry like 
time reversal or simplex, the block structure of the ma- 
trices involved fixes the sign. This has been discussed, 
for instance, in some recent applications of angular mo- 
mentum projection (AMP) using axially symmetric and 
time reversal preserving intrinsic wave functions 0, [J]. 
To move forward, HFB wave functions that do not have 
any spatial symmetry (triaxial) and also breaking time 
reversal symmetry have to be considered in order to in- 
corporate K ^ configurations. This is the case to de- 
scribe, for instance, the ground state of odd-A nuclei. 
For the usual time reversal breaking (TRB) mean field 
wave functions, the simplex symmetry endows the HFB 
amplitudes U and V with a common bipartite structure 
and the usual arguments used to extract out the sign 
of the overlap apply. However, when full triaxial an- 
gular momentum projection of HFB intrinsic states 0] 
is considered, the simplex symmetry is no longer pre- 
served in the evaluation of rotated overlaps and the de- 
termination of the sign becomes more difficult. A gen- 
eral solution to the sign problem was given in Ref. 0], 
where it was shown that the overlap, including the sign, 
can be computed from the pairwise degenerate eigen- 
values of a non-hermitian matrix. Handling the eigen- 
values of non-hermitian matrices is a difficult task 0], 
that increases its complexity if the pairwise degenerate 
eigenvalues have to be obtained numerically without any 
symmetry enforcing degeneracy, as is the case with HFB 
wave functions breaking simplex. Neergard's method has 
been used along with small configuration spaces [8| but 
in the majority of the calculations continuity arguments 
are used (see Refs. 0,0, E3] for recent examples) in spite 
of the difficulties with that procedure. The same sign 
problem is also present in the evaluation of the trace of 



statistical density operators [Til, E2]- In this case, how- 
ever, Neergard's method has not been implemented up to 
date, leaving as the only choice the continuity method in 
such finite temperature calculations. The same difficulty 
also applies to the recently proposed method to compute 
multiquasiparticle overlaps that relays on the statistical 
Wick's theorem [l3|]. Recently, the group structure 
of the unitary Bogoliubov transformation has been dis- 
cussed, as well as its implications in the relative phase 
between two HFB wave functions. However, its practical 
implications are still unclear. 

In this paper, I will introduce a new way to compute 
the overlap of two HFB wave functions based on the con- 
cept of fermion coherent states [H, EE, G3] ■ The new for- 
mula involves a quantity similar to the determinant called 
pfaffian of a skew-symmetric matrix. The advantage of 
the proposed method is that the numerical evaluation of 
the pfaffian is simple and lacks the problems previously 
mentioned about pairwise degenerate eigenvalues. An- 
other advantage of the present formulation is its applica- 
bility to the evaluation of the trace of density operators 
like the ones found in applications of the auxiliary-field 
shell model Monte Carlo[ll| or symmetry restoration at 
finite temperature[T3|. A reliable determination of the 
sign of the norm can also be useful in to order to pin 
down the location of the zeros of the HFB overlaps 
This determination would eventually be useful to get rid 
of the so called "pole problem" that plagues present be- 
yond mean field calculations. 

II. OVERLAPS AND TRACES 

A. Preliminaries 

Let \4>q) and \4>\) be two HFB wave functions defined 
in terms of a set of single particle creation and anni- 
hilation operators and at that are assumed to be re- 
lated by hermitian conjugation and also to satisfy fermion 
commutation relations. The HFB wave functions, in the 
Thouless representation 0, HE], are given by 

i^) = ex puE M ^^4)io) a) 

V kk> / 



where the skew-symmetric matrices 



2 



are defined in terms of the f/j and Vi coefficients of the 
Bogoliubov transformations defining the HFB wave func- 
tions and |0) is the true vacuum. The arbitrary phase 
that can always be associated with a vector state in quan- 
tum mechanics has been implicitly fixed in the definition 
of Eq. (HJ by requiring (0\4>i) = 1. Ways to enforce 
this normalization for general HFB wave functions are 
discussed, for instance in Refs. 0, In the event 

of having (0\4>i) — (as a consequence of divergent 
and/or zero occupancies) the best practical strategy is to 
use another reference wave function instead of the true 
vacuum |0). The new reference HFB wave function \4>) 
has to be conveniently chosen as to stay close to both 
\<j)i) (for instance by taking a wave function with similar 
deformation parameters as those of \4>i)). The matrices 
MW referred to \4>) can be straightforwardly computed 
in terms of the previous quantities and the Bogoliubov 
transformation amplitudes of the reference state. In the 
rare event of not finding a convenient reference wave func- 
tion |</>) a regularization procedure to handle the diver- 
gent AfW matrix elements (or the zero occupancies) is 
in order. In this case, the expressions get more involved 
and a detailed account is deferred to a forthcoming pub- 
lication. Another way to deal with that problem is pre- 
sented in Ref. [f| but the resulting expressions are rather 
involved. 

Let me now introduce fermion coherent states |z), 
which are parametrized in terms of the anticommutingel- 
ements Zk and zl of a Grassmann algebra [l5l.fl^.fl7l.ll8l| 
and fulfilling the equations 



B. Evaluation of the overlap 



To compute the overlap (4>o | 
Eq. (@]) is inserted to obtain 



l), the closure relation of 
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Using now Eqs. ([2j) and one arrives to 



(5) 

where the property |(0|z)| 2 = 1 is used. The integral is 
of the Gaussian type but for Grassmann variables. The 
techniques to evaluate this kind of integrals can be found 
in many textbooks flR 16, [l7| but its evaluation will be 
carried out explicitly here. The reason is that in order 
to determine the sign of the norm we have to be careful 
with some intermediate steps. The above integral can 
be written in a more compact way by introducing the 
bipartite skew-symmetric matrix 



k'k 



and the vector of Grassmann variables z M = {zy, Zk>) as 



Wi) = [H(dz* k dz k )ei^>z^>^ 



(6) 



a k \z) = z k \z) 



(2) 



and 



(3) 



From the above definition is clear that |z) is a right eigen- 
state of the annihilation operator ak with eigenvalue z k 
whereas (z| is a left eigenvector of at with eigenvalue 
zl (The notation used for the members of the Grassmann 
algebra is the usual one but can be a little misleading as 
zl is not connected to Zk by complex conjugation). The 
coherent states satisfy a closure relation 



d/i(z)|z)(z| 



(4) 



where the metric of the integral is given by dfi(z) = 
e ~ z z life dz* k dzk. These and other relevant definitions 
and properties of fermion coherent states can be found 
in many textbook or in the original literature (TEI . 



The skew-symmetric matrix M can always be trans- 
formed [iji] to canonical form by means of a unitary 
transformation U 



= U 



( o 



-Pi 



V o 



pi \ 












-Pn 









Pn 
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U T = UM C U T 





and the Pi,.. . ,Pn coefficients of the "canonical form" 
of the matrix M are real and positive. Introducing now 
the new Grassmann variables r]^ = Y1^(^ T )^^ Z ^ the 
exponent in the integrand of Eq. |(5]) becomes 



2 ^ 

mm' 



N 

E 

k=l 



PkVkVk 



which is straightforward to integrate. The Jacobian of 
the transformation can be shown to be simply det(U T ) — 
det(U). The remaining integrals can be performed easily 
being the result 



drfdrje 



0V*v — 



-P- 
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The final expression for the overlap is then 



N 



6o|0i> = (-If det(U)Y[(3 k . 



k=l 



This expression can be cast in terms of the pfaffian 
of a skew-symmetric matrix. The pfaffian of a skew- 
symmetric matrix (see, for instance, [20|]) is a number ob- 
tained out of the matrix elements of the skew-symmetric 
matrix in a way quite similar to the one used to define 
the determinant (see Appendix A for details and prop- 
erties used below). The connection between the product 
of /3i's and the pfaffian is a consequence of Eq. (| A2[) 
and reads Uk=i Pk = (~l) N( - N ' 1)/2 pf{M c ) where pf(M c ) 
obviously denotes the pfaffian of M c . Using the prop- 
erty GQJ pf(M) = pf(UM c U T ) = det(L/)pf(M c ) the final 
result is obtained, 



(00 101 ) = SArpf(M) = S N pf 



MW -1 
1 -M(°)' 



(7) 



where sn = (— 1) N ( N+1 '' 2 , To make the connection 
with the standard formula for the overlap Q the rela- 
tion pf(A) 2 = det A is used (and this is here where the 
sign is lost) to write 







det 



i 



-M(°)' 



1/2 



(8) 



This expression reduces, by using standard formulas for 
the determinant of a bipartite matrix (see below) , to 



6 o |0i) - (det(l -M(°)«mW)) 



1/2 



(9) 



which is the usual expression for the norm (Onishi for- 
mula). Please notice that in going from Eq. (7]) to Eq. 
([8]) the sign present in the first equation is lost as a con- 
sequence of the writing of the square of the pfaffian as 
a determinant. Also signs appearing in the manipula- 
tions needed to obtain Eq. |(9]) have been neglected. We 
clearly see that the sign problem appears in the standard 
formulas because of the wrong implicit use of the above 
relation between the pfaffian and the determinant. 

If both HFB wave functions |0o) and |0i) share a com- 
mon discrete symmetry like simplex or time reversal, then 
the matrices defining them can acquire a common 
block structure 







M 
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-M WT 



that can be used to simplify the result of Eq. (0. In this 

case 
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By exchanging blocks 2 and 4 we obtain 



pf(M) = (-ir pf 
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to give 







6o|0i)=det( i, M ' ) =det(l + Af ( ° )+ M (1) ). 
\M 1 J 

(10) 



C. Evaluation of statistical traces 

Now I turn to the evaluation of the trace of density 
operators. In the^ statistical HFB theory the statistical 
density operatorPis given by the exponential of a one 
body operator V = exp[| ^f^TZ^u] where 7 M is a 
shorthand notation for {fix, . . , ,0n, (iy , . . . , (3^) and 1Z 
is a skew-symmetric matrix of dimension 2N character- 
izing the density operator (see Ref [13] for details). An- 
other way to characterize the density operator is to de- 
fine how it transforms quasiparticle creation and annihi- 
lation operators V~ lr y^V = Ylv^f 1 "^" where the matrix 
T = exTp(alZ) and cr^ v — {7^,71,}. Introducing the bi- 
partite structure of T 



T : 



T21 T22 



the Balian and Brezin's decomposition [2l| of V is given 

by 



V 



with X — T12T22 and Z — T 22 T21 skew-symmetric (as 
a consequence of the relation T T aT — a that T satisfies) 
and exp(— Y) = T^. To evaluate the trace of T> using 
fermion coherent states we have to use the formula [ItJ 

Tt(D) = J dfi(z){-z\V\z) (12) 

where |z) are again a set of fermion coherent states but 
chosen this time as eigenstates of the quasiparticle anni- 
hilation operators i.e. (3i\z) — Zi\z). Using Eq. ifTTj) 
the evaluation of the overlap between the fermion coher- 
ent states gives 

(-z\V\z) = e-^^e^n z * Xi ^e^ ZiZ ^ 

To evaluate the remaining overlap the standard re- 
sult exp (j^ij PtYijPj) \z) = \e Y z) used together with 
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-z|z') = exp(-z*z') (see Refs. [H, [H, ) gives 



Combining this result with Eq. (|T2[) . the following inte- 
gral is obtained 

Tr(T>) = e-i^W J JJ (dz* k dz k ) 2w V m m'm** 

k 

where the same notation as in Eq |(6]) is used. In this case 

X -{e Y + t)' 
(e Y + l) T ' Z 

Applying the same considerations as in the evaluation of 
the overlap we finally arrive to 

Tr(D) = s N exp f~Tr[Y]j pf(M) 

where — (—l) N ^ N+1 ^/ 2 . Taking into account the re- 
lationship between X, Z and Y and the blocks of the 
matrix T the above result can be expressed as 



Tr(£>) = sjv(detT 22 ) 
x pf 



1/2 



T'l2'7 1 2 2 



22) 
T21 



The introduction of (detT 22 



,V2 



in place of 



exp (— iTr[Y]) can lead to the (right) conclusion 
that a sign indeterminacy has been introduced in the 
expression of the trace. The definition of V in terms 
of the transformation matrix T leaves a phase open 
in the definition of the density operator which is also 
present in the expression of Eq. fTTj) . A way to fix the 
phase is to require some condition like, for instance, the 

A 1/2 

realness and positiveness of (0o|f| ( / , o) = (detT 22 ) ' 
where \<po) is the vacuum of the quasiparticle operators 

Pi entering in the definition of V. This condition implies 

1/2 

the replacement of (detT 22 ) by its modulus. Using 
property l|Al|) of the pfaffian the final result is obtained 



Tr(2>) 



|detT 2 



[1/2 



S N - 



x pf 



det T 22 
T12T22 
(T22 + t) 



L 22 T 1) 
T21T22 



(13) 
(14) 



This result is apparently quite different from the standard 
one of [lll.[l3|. but after some tedious manipulations (see 
appendix B) one can obtain the usual result. 



III. CONCLUSIONS 



HFB wave functions. The result given in terms of pfaffi- 
ans is simpler to implement than previous considerations 
[6] based on pairwise degenerate eigenvalues of a general 
matrix and it is free from the uncertainties of other meth- 
ods based on continuity arguments. Indications on how 
to evaluate efficiently the pfaffian are also given. Hope- 
fully, this new method will help to simplify the implemen- 
tation of ambitious projects like triaxial angular momen- 
tum projection. On the other hand, the method used is 
straightforwardly extended to the evaluation of the sign 
of the trace of statistical density operators which is a new 
result not considered previously in the literature. 
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Appendix A: DEFINITION, BASIC PROPERTIES 
AND NUMERICAL EVALUATION OF THE 
PFAFFIAN 

The pfaffian of a skew-symmetric matrix R of dimen- 
sion 2N and with matrix elements is defined as [2Gj| 



pf(i?) 



1 1 

2" n! 



1-2 ' «3U ' «5*6 



■ ?"2n-l,2n 



Perm 



where the sum extends to all possible permutations of 
ii, . . . , i 2n and e(P) is the parity of the permutation. 
For matrices of odd dimension the pfaffian is by def- 
inition equal to zero. As an example, the pfaffian of 
a 2 x 2 matrix R is pf(i?) = ri 2 and for a 4 x 4 one 
pf(i?) = ri 2 r 3 4 — ri 3 r 2 4 + 7*14^3. Similarly to the case 
of determinants, exchanging rows i and j and the same 
time columns i and j, multiplies the pfaffian by minus 
one. Other useful properties of the pfaffian are 



pf(P T RP) = det(P)pf(-R), 



.. , R 

pf l -R T 



) =(-l)^ A, - 1 )/ 2 det( J R) 



pf (R) = (det R) 



1/2 



(Al) 



(A2) 



(A3) 



A useful formula to compute pfaffians of small or simple 
matrices is 

pf(iZ) = E(- 1 ) l+J " 1 ^P f ( i? y) (A4) 



I have used the technique of fermion coherent states 
to compute unambiguously the sign of the overlap of two 



where Rij is the pfaffian-minor obtained by eliminating 
from R the two rows and two columns i and j. 



5 



The pfaffian of a complex skew-symmetric matrix 
R of dimension 2N is evaluated numerically by first 
reducing the matrix to tridiagonal form Rt- This 
reduction is accomplish by means of a set of 2(N- 
1) successive Householder transformations Pi exactly 
in the same way as in the standard reduction of a 
symmetric matrix to tridiagonal form pj. We have 
P 2 {N-1) ■ ■ ■ P 2 PiRP?P? ■ ■ ■ PrffN-D = Rt with 



/ n 

-n 



Rt 





r 2 



-r 2 












\ 



r 2 N-i 
-r 2N -i / 



where the special structure of a skew-symmetric and 
tridiagonal matrix is evident. Using now Eq. 1|A1[) we ob- 
tain det(Pi) . . .det(F 2 (A r -i))pf(^) = pf(-RT). As House- 
holder matrices are hermitian, unitary and have determi- 
nant det(Pi) = — 1 we finally arrive to pf(-R) = pf(i?r)- 
To evaluate the pfaffian of the tridiagonal matrix we use 
the minor expansion of Eq. (| A4|) that gives pf (Rt) 

t-tN 

nr 3 . ..r 2 N-i 



lL=l r 2i-l. 



Appendix B: DERIVATION OF THE STANDARD 
FORMULA FOR THE TRACE 



Using Eq. (| A3[) the pfaffian of M is written as the square 

/ -N.l/2 

root of its determinant, pfM = (det Ml . By exchang- 
ing rows and columns conveniently 



detl 



(-) 



\ 1 (T 22 + 1) T 21 T£ 2 
T 12 T£ -(T 2 t 2 + U) 



and applying the formula of a bipartite determinant 

det ( C D ) = dct A Act ( D ~ CA ^ B ) the following re- 
sult is obtained 

dctM = detT 22 det(l + T 22 ) 

x det (t + (T 2 T 2 ) _1 + T 12 T^(l + T 22 )- l T 2 , 

But Tn = (T^)' 1 + T 12 T^T 21 and 
T22 1 [(1 + T 22 )- x - i] = -(1 + T 22 y x so that 

detM = detT 22 det(l + T 22 ) 

x dct (1 + Tn - T 12 (t + T 22 y x T 21 ) 

' (Tn + 1) Ti 2 
T 2i (T 22 + l) 
When the pfaffian is written as the square of the deter- 
minant the sign is lost and therefore phases are irrelevant 
in the derivation. Taking all this into account the result 



det T 22 det 



In this appendix the standard result of [12J for the trace 
of a density operator is deduced from Eq. lfl3|). I start 
considering 



P\-iP 22 



22 



-(7f 2 + l) 
T2iT 22 



Tr(T>) 



det 



(Tn + 1) Ti 2 
T21 (T 22 + 1) 



1/2 



is obtained up to a sign, which is the sought formula of 
Reffll. 
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